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LONG-TERM  GOALS 

Development  of  Regional  Oceanic  Modeling  System  (ROMS)  with  emphasis  on  the  ability  to 
simulate  realistic,  highly-turbulent  flows  in  both  basin-  and  local  regional  scales,  at  resolutions 
sufficient  to  capture  submesoscale  phenomena;  analysis  and  understanding  the  underlying  physical 
processes;  improving  parameterizations  of  unresolved  processes,  applicable  numerical  algorithms, 
relevant  computer-science  techniques  needed  for  efficient  parallel  code  adapted  for  modem 
computing  environment. 

OBJECTIVES 

The  scientific  question  is  how  the  eddies  control  the  currents  by  the  eddy-induced  momentum  and 
buoyancy  fluxes,  turbulent  mixing,  and  their  bottom  form  stress  (pressure  force)  and  bottom  boundary 
layers  -  all  the  aspects  associated  with  turbulent  flows  over  steep  topography  in  the  presence  of 
stratification.  This  is  done  in  the  context  of  western  boundary  currents,  the  Gulf  Stream  and  the 
Kuroshio,  as  well  as  the  western  Pacific  (Solomon  Sea).  The  intend  is  obtain  realistically  accurate 
separation  and  the  subsequent  path  of  the  Gulf  Stream,  the  Kuroshio  path  from  Taiwan  to  Japan  with 
subsequent  separation,  as  well  as,  provide  analysis  of  modeling  sensitivities  to  input  parameters  and 
forcing,  and,  wherever  possible,  to  understand  and  give  a  physical  explanation  of  the  underlying 
mechanisms.  This  involves  validation  against  real-world  data  in  cooperation  with  William  S.  Kessler 
and  Hristina  Hristova  from  PMEL  (Solomon  Sea),  and  Satoshi  Miarai,  Taichi  Sakagami  from 
Okinawa  Institute  of  Science  and  Technology  (Kuroshio). 

APPROACH 

The  primary  method  is  to  obtain  numerical  solution  on  sets  of  progressively  refined  grids,  starting 
with  basin-wide  eddy  permitting  resolutions  (although  substantially  finer  than  that  used  in  climate 
modeling),  and  downscaling  it  to  Ax  ~  1  km  or  less  which  enables  adequate  resolution  of  mesoscale 
and  submesoscale  eddies  generated  by  instability  and  bottom  topography  effects.  This  is 
accompanied  by  development  of  modeling  codes,  numerical  methods,  and  the  associated  pre-  and 
post-processing  infrastructure,  analysis,  visualization  that  allows  prioritizing  the  tasks  directly 
necessary  for  achieving  the  objectives  stated  above,  but,  at  the  same  time,  traditionally  maintaining 
the  3-way  balance  among  physical  applications  (including  introducing  new  submodels),  numerical 
methods,  and  computer  science  that  has  been  established  throughout  ROMS  development  history. 
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Additionally,  we  design  idealized,  process-study-oriented  simulations  with  the  intent  of 
understanding  of  physical  phenomena,  and  as  control  tests  for  algorithm  verification. 

WORK  COMPLETED 

During  the  last  year  we  worked  on  decadal  Pacific  and  Atlantic  circulations,  the  Kuroshio,  and  the 
Gulf  Stream;  mesoscale  eddy  buoyancy  fluxes;  submesoscale  surface  fronts,  filaments,  and  eddies; 
topographic  current  separation,  form  stress,  and  submesoscale  vortex  generation;  Our  work  on 
isoneutral  diffusion  for  tracers  is  now  published,  Lemarie,  et  al. (201 2a, b),  however  most  recently  we 
have  revisited  the  relevant  parts  of  the  code  with  the  intent  of  reducing  its  computational  cost,  and  for 
better  integration  with  the  rest  of  the  model.  On  the  new  algorithmic  direction  we  implemented 
adaptively  implicit  vertical  advection  to  alleviate  the  associated  CFL  limitations,  which  has  been  a 
growing  impediment  for  high-resolution  simulations,  as  well  as  for  modeling  tides.  The  code  now 
utilizes  dual  (shared-  and  distributed-memory)  parallelization  (implemented  via  Open  MP  and  MPI, 
respectively,  with  an  option  to  use  both  at  the  same  time)  based  on  subdomain  decomposition  for  both 
with  two-level  hierarchy.  We  also  made  an  effort  to  parallelize  pre-processing  software  tools 
wherever  justified  by  their  computational  cost.  A  substantial  effort  on  analysis  and  visualization 
techniques  was  performed  by  Jonathan  Gula  that  involves  transition  from  Matlab  to  Python-based 
software  and  graphics. 

RESULTS 

We  present  a  few  highlights  for  this  project.  The  publications  list  (papers  from  2012  up  to  ones  likely 
to  be  submitted  in  2013)  provides  a  view  of  the  finalized  results  across  all  our  ONR  projects. 

Eddy  buoyancy  fluxes:  It  is  now  widely  understood  that  mesoscale  eddy  fluxes  play  a  central  role  in 
the  oceanic  general  circulation,  requiring  that  they  either  be  parameterized  (as  in  most  climate 
simulations)  or  adequately  resolved  with  grid  scales  of  Ax  =  10  km  or  smaller.  Nevertheless,  their 
spatial  heterogeneity  and  long  time  scales  of  dynamical  equilibration  make  them  difficult  to  measure 
accurately  and  challenging  to  interpret  in  models.  A  dynamical  interpretation  is  made  of  the 
mesoscale  eddy  buoyancy  fluxes  in  the  Eastern  Boundary  Currents  of  California  and  Peru-Chile, 
based  on  regional  equilibrium  ROMS  simulations.  The  eddy  fluxes  are  primarily  shoreward  and 
upward  across  a  swath  several  hundred  km  wide  in  the  upper  ocean  (Fig.  1);  as  such  they  serve  to 
balance  mean  air- sea  heating  and  cooling  by  coastal  upwelling.  In  the  stratified  interior  the  eddy 
fluxes  are  consistent  with  the  adiabatic/isopypcnal  hypothesis  associated  with  a  mean,  eddy-induced 
velocity  advecting  mean  buoyancy  and  tracers  in  a  materially  conservative  way  that  extracts  potential 
energy  from  the  mean  flow.  Furthermore,  with  a  suitable  gauge  choice,  the  horizontal  fluxes  are 
almost  entirely  aligned  with  the  mean  horizontal  buoyancy  gradient,  consistent  with  the  advective 
parameterization  scheme  of  Gent  and  McWilliams  ( JPO ,  1990).  The  associated  diffusivity  is 
surface-intensified,  matching  the  vertical  stratification  profile  below  the  boundary  layer.  The  fluxes 
span  the  across-shore  band  of  high  eddy  energy  (hundreds  of  km  wide),  but  their  along-shore 
structure  is  unresolved  due  to  statistical  sampling  limitations.  In  the  surface  layer  the  eddy  flux  is 
significantly  diabatic  with  a  shallow  eddy-induced  (Lagrangian)  circulation  cell  and  down-gradient 
lateral  diapycnal  flux.  These  effects  are  summarized  in  Fig.  2.  The  dominant  eddy  generation  process 
is  baroclinic  instability,  but  there  are  significant  regional  differences  between  the  different  upwelling 
systems  in  the  flux  pattern  and  the  associated  k  profiles  that  are  not  consistent  with  simple  baroclinic 
instability  theory  (Colas  et  al.,  2013a). 
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Dynamical  balance  of  the  Gulf  Stream:  Identifying  the  dynamics  responsible  for  the  Gulf  Stream 
separation  and  control  of  its  path  has  been  a  long-standing  challenge  in  oceanography.  The  dynamical 
balance  of  the  Gulf  Stream  is  investigated  using  very  high-resolution  ROMS  simulations  highlighting 
the  important  role  played  by  eddy-driven  flows  and  the  interactions  with  complex  topography. 

The  bottom  pressure  torque,  which  derive  from  the  twisting  of  the  force  that  the  bottom 
topography  exert  on  the  ocean,  is  the  main  ingredient  enabling  the  return  flow  of  the  wind-driven 
transport  in  western  boundary  currents.  The  bottom  pressure  torque  strongly  reflects  the  regional 
topography  as  seen  on  Fig.  3.  The  planetary  vorticity  term,  related  to  meridional  transport  of  fluid,  is 
balanced  by  the  sum  of  the  bottom  pressure  torque  and  the  non-linear  advective  term  in  the  barotropic 
vorticity  balance  equation.  The  non-linear  advective  term  redistributes  vorticity  locally  but  spatially 
averages  to  zero  for  areas  larger  than  «  2°  x  2°.  Hence  the  bottom  pressure  torque  is  the  term 
providing  the  overall  positive  input  of  vorticity,  balancing  the  negative  input  by  anticyclonic  wind 
curl  on  the  scale  of  the  subtropical  gyre.  Differences  between  the  two  bottom  panels  of  Fig.  3, 
especially  in  the  Charleston  Bump  region,  are  due  to  the  contribution  of  the  bottom  stress  curl.  Fig.  4 
shows  how  the  Gulf  Stream  path  is  directly  linked  to  the  Bottom  Pressure  Torque,  and  its  related 
eddy-driven  abyssal  circulation,  through  vertical  fluxes  of  vorticity. 

Energy  conversions  at  the  topographic  features,  especially  the  Charleston  Bump  and  at  the 
separation  point  after  Cape  Hatteras,  show  that  they  are  triggers  for  instability  processes  and 
baroclinic-barotropic  energy  conversions.  The  eddy  potential  to  eddy  kinetic  energy  ( w'b ' ,  Fig.  5) 
shows  a  strong  baroclinic  conversion  right  after  separation  and  near  the  Charleston  Bump,  with  very 
small  or  negative  signals  in-between.  The  mean  kinetic  to  eddy  kinetic  energy  due  to  the  horizontal 
and  vertical  Reynolds  stresses  ( KmKe  =  HRS  +  VRS,  Fig.  5)  also  shows  a  barotropic  conversion  in 
the  same  places.  The  topographic  constraint  has  a  stabilizing  effect  on  the  flow  upstream  from  the 
Charleston  Bump.  This  constraint  is  reduced  as  the  stream  is  deflected  seaward  by  the  Bump  where 
the  topography  is  deeper  and  a  relatively  large  conversion  of  mean  energy  to  eddy  energy  can  happen. 
The  topographic  constaint  is  renewed  between  the  Bump  and  Cape  Hatteras  while  the  Gulf  Stream 
flows  northeastward  and  requires  that  some  eddy  energy  be  converted  back  to  mean  energy.  The 
strongest  eddy  generation  happens  after  separation  of  the  Gulf  Stream  due  to  the  increased  transport 
and  the  confluence  of  the  flow  which  acts  to  sharpen  the  cross  frontal  gradients  of  the  mean  flow. 

Adaptively-implicit  vertical  advection  for  momentum  and  tracers:  As  any  oceanic  model  with  an 
Eulerian  vertical  coordinate,  ROMS  is  subject  to  the  CFL  limitation  associated  with  vertical 
advection.  When  the  goal  is  to  produce  a  simulation  with  horizontal  resolution  less  than  a  few  km 
this  limitation  gradually  becomes  the  most  restrictive  one  (depending  on  the  topographic  features 
and/or  whether  or  not  the  tides  are  part  of  the  simulation).  For  sub- 1  km  grids  it  becomes  so  dominant 
that  the  imposed  time  step  limitation  is  several  times  smaller  than  that  due  to  horizontal  advection 
and/or  internal  wave  phase  speeds.  Detailed  investigation  of  the  associated  “hot  spots”  (characteristic 
locations  on  the  model  grid  where  numerical  instability  of  the  explicit  code  occurs  first )  reveals  large 
vertical  CFL  always  occur  near  topographic  features,  where  buoyancy  stratification  is  weak  or 
vanishing,  but  not  necessarily  in  the  shallowest  areas  where  vertical  grid  spacing  is  the  smallest  due 
to  topography-following  coordinate. 

While  implicit  advection  schemes  offer  a  relief  from  CFL  limitation,  their  drawbacks  are  well 
known:  unavoidable  and  potentially  large  dispersive  errors  increasing  with  CFL,  and  depending  on 
the  detail  of  time  and  space  discretization,  large  numerical  viscosity  as  well  ( cf Shchepetkin  & 
McWilliams,  2008).  In  contrast,  the  explicit  vertical  advection  schemes  of  ROMS  are  designed  to  be 
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high  order  in  space  (4th-centered  or  compact  based  parabolic  spline  fits),  which  makes  it  not  feasible 
to  design  an  implicit  version  of  comparable  accuracy. 

To  overcome  the  dilemma  we  pursue  an  adaptive  approach  where  the  scheme  remains  explicit  (as 
in  the  original  code)  everywhere  except  where/when  local  vertical  velocities  exceed  a  threshold  close 
to  (but  below)  the  explicit  stability  limit.  Once  this  happens,  a  gradual  transition  toward  an  implicit 
scheme  begins  via  Courant-number-dependent  weighting  algorithm.  As  we  are  not  aware  of  any 
prototype  of  such  approach  published  in  the  literature,  we  present  a  few  detail  here.  Vertical 
advection  fluxes  for  the  tracer  or  velocity  fields  are  discretized  in  such  a  way  that  their  computation 
involve  the  advected  field  at  the  new  time  step,  q^+1  which  is  yet  unknown, 

FCt+^  =  Wk+^-£{ql+\qliUT1'\ql?f,.-)  ■  (1) 

This  can  be  rearranged  by  splitting  vertical  velocity  into  two  parts, 

W4+1/2  =  w£>1/a  +  Vt  — 0,1 . AT  (2) 


where  W^y2  is  used  to  compute  terms  involving  q%+1,  q%±l  only  (/.<?.,  implicit  part),  while  W^y2  for 

the  remaining  ++  +  (h±/'2 >  ■■■  (half-integer  time  index  n  +  l/‘>  means  that  these  are  provisional 
values  from  ROMS  predictor-corrector  stepping,  Sec.  4  in  Shchepetkin  &  McWilliams,  2005).  Then, 
-terms  are  computed  via  standard  algorithm  within  the  r.h.s.,  while  -terms  are  integrated 
into  the  vertically  implicit  operator.  Assuming  upstream  treatment  of  the  implicit  part, 
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the  combined  implicit  advection-diffusion  system  becomes 
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where  the  prime  in  vhsk  means  that  the  usual  r.h.s.  computed  by  ROMS  code  for  the  corresponding 
equations,  except  the  replacement  W4+i/2  — >  Ak+ y2  is  vertical  viscosity/diffusion  coefficient 

(including  the  stabilization  terms,  Sec.  3  from  Lemarie,  et  al.( 2012b),  in  the  case  when  isoneutral 
lateral  diffusion  is  used).  The  above  system  takes  into  account  kinematic  b.c.  at  surface  and  bottom, 
Wn± i/2  =  kki/2  =  0,  bottom  no-flux  b.c.  for  tracers  (Eq.  (6)  for  the  momentum  equation  contains  an 
extra  term  associated  with  bottom  drag,  which  also  treated  implicitly).  As  in  the  original  ROMS  code, 
the  new  algorithm  has  simultaneous  conservation  and  constancy  preservation  properties  for  tracers, 
despite  the  fact  that  grid  box  heights  change  due  to  changing  free  surface,  Hk+1  4  Hk. 

The  splitting  in  (2)  works  as  follows:  At  first  the  vertical  velocity 
114+  '/a  is  computed  the  standard  way,  and,  in  addition  to  that,  also 
computed  is  the  finite-volume  Courant  number  Cuk  (defined  at  every 
grid  box  Hk  as  the  sum  of  fluxes  outgoing  from  Hk,  normalized  by 
the  time  step  size  At  and  grid-box  volume).  The  explicit  part  is  then 
computed  as 


w(e)  Wk+ 1/2 

k+ 1/2  f(C*)  ’ 


C*  =  Cuk  if  114+1/2  >  0 

C*  =  Cuk+i  if  114+1/2  <  0 


while  the  implicit  is  the  remainder, 


=  114+1/2  -  w,(e) 


fc+1/2 


fc+1/2 


The  limiting  function  /  =  f(Ou )  is  shown  on  the  left,  and  consists 
of  three  segments:  constant,  parabolic,  and  linear,  smoothly  matched 
to  each  other.  The  two  parameters  Oum in  and  Cumax  control  the 
threshold  below  which  the  algorithm  is  fully  explicit  as  is  the 
original  ROMS  code,  and  the  maximum  allowed  Courant  number 
for  the  explicit  part.  The  practical  values  we  select  in  our 
computations  are  Cum in  =  0.2  and  Oumax  =  0.5,  which  is  based  on 
consideration  of  numerical  stability,  and  accuracy  as  well. 

The  motivation  for  using  upstream  discretization  for  the  implicit  part  comes  from  the  fact  that  while  it 
is  diffusive,  it  is  monotonic  and  has  smaller  numerical  dispersion  in  its  truncation  term  than  the 
centered  scheme  (due  to  avoidance  of  spatial  averaging).  This  leads  to  a  well  posed,  diagonally 
dominant  discrete  system.  This  choice  is  further  justified  by  the  observation  that  in  practical  model 
solutions  large  vertical  velocities  occur  only  in  places  with  vanishing  (or  even  unstable)  stratification 
and,  consequently,  already  large  mixing  set  by  the  vertical  parameterization  scheme. 


To  demonstrate  the  viability  of  the  proposed  algorithm  we  present  test  results  from  gravitational 
adjustment  “lock-exchange”  problem,  Fig.  6.  This  problem  is  known  for  generating  sharp  fronts  with 
the  resulting  vertical  velocities  playing  the  dominant  factor  in  CFL  stability.  The  adaptively  implicit 
algorithm  allows  dramatic  increase  of  the  allowed  time  step.  For  the  smallest  settings  of  At,  adaptive 
and  explicit  solutions  are  identical  (which  is  expected),  while  for  the  largest  At  adaptive  solution 
becomes  more  similar  to  fully  implicit  (backward  Euler  in  time,  upstream  in  space).  Also  note  the 
progressive  delay  in  the  front  propagation  for  the  largest  At  -  neither  adaptive,  nor  fully  implicit 
scheme  is  expected  to  be  accurate  at  this  regime  (At  =  240. ..510s),  but  still  adaptive  shows  slightly 
less  delay  and  less  mixing. 
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TRANSITIONS 


ROMS  modeling  code  and  the  associated  supporting  software  are  generally  available  as  open  source. 
We  also  made  available  some  of  our  solutions  for  the  analysis  by  a  third  party. 

RELATED  PROJECTS 

Wave-current  interaction  (in  cooperation  with  Yusuke  Uchiyama,  University  of  Kobe,  Japan); 
Biogeochemistry  modeling  (in  cooperation  with  Curtis  Deutsch,  Jun-Hong  Liang,  and  Hartmut 
Frenzel,  from  University  of  Washington) 
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Figure  1:  Annual  and  along-shore  averages  of  across-shore  (left)  and  vertical  (right)  eddy  buoyancy 
flux  [in2  s~3]  for  the  Peru  (top),  Chile  (middle),  and  California  (bottom)  current  systems.  White 
contours  are  the  mean  buoyancy  field,  and  the  gray  dashed  line  is  the  boundary-layer  depth  diag¬ 
nosed  by  the  K-Profile  Parameterization.  Black  contours  are  vertical  flux  w'b'  =  0.  The  eddy  flux 
patterns  are  generally  shoreward  in  u'b'  and  upward  in  w'b',  acting  to  balance  the  air-sea  heating 
and  mean  upwelling  circulations.  While  there  is  general  qualitative  similarity  across  the  regions, 
there  are  noticeable  regional  differences  in  stratification  and  fluxes.  (Colas  et  al.,  2013a) 
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Figure  2:  3D  Schematic  representation  of  the  eddy  effects  on  the  mean  buoyancy  field  decomposed 
between  adiabatic  eddy-induced  advection  and  diabatic/diapycnal  flux.  ( The  geometry  is  for  Peru 
with  its  Andes  mountains.)  The  mean  stratification  is  indicated  by  the  contours  in  the  foreground 
vertical  plane  as  well  as  by  the  colors.  In  the  surface  boundary  mixed  layer,  the  adiabatic  advective 
component  is  associated  with  the  restratication  tendency  of  submesoscale  fronts  with  a  lateral  scale 
0(1  km).  It  is  represented  by  its  associated  streamfunction  ( black  closed  contours  in  the  foreground 
vertical  plane).  The  diabatic  component  acts  to  smooth  out  surface  buoyancy  extrema  and  is  shown 
as  sinuous  arrows  in  the  top  plane.  Interior  diabatic  fluxes  represented  by  white  arrows  in  the 
foreground  vertical  plane  are  found  to  be  negligible  in  our  mesoscale-resolving  solutions.  Interior 
adiabatic  advective  fluxes  associated  with  mesoscale  eddies  ( whose  associated  streamfunction  is 
shown  with  gray  closed  contours  in  the  foreground  vertical  plane)  oppose  Ekman-induced  transport. 
Both  adiabatic  components  can  be  seen  as  the  manifestation  of  Available  Potential  Energy  release 
by  baroclinic  instability  and  frontogenesis  near  the  surface.  (Colas  et  al.,  2103b) 
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Figure  3:  Time-mean  of  (a)  Bottom  Pressure  Torque  J^b0’h\  (b)  horizontally  smoothed  Bottom  Pres¬ 
sure  Torque,  (c)  sum  of  the  Bottom  Pressure  Torque  and  the  non-linear  advective  terms  in  the 
barotropic  vorticity  balance  equation  —  A  and  (d)  the  planetary  vorticity  term  — V.(fu)  for 

18  years  of  a  2.5  km  resolution  ROMS  simulation.  Units  are  m.s~2.  Topography  is  shown  in  black 
lines  and  mean  barotropic  streamfunction  in  green  lines  (contours  every  10  Sv). 
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Charleston  Bump 


Figure  4:  Mean  vertical  advection  of  absolute  vorticity  for  a  18-year  2.5 km  resolution  ROMS 
simulation.  Blue  contour  is  —0.5 e~8ms~2  (downward)  and  red  contour  is  +0.5 e~8ms~2  (up¬ 
ward).  The  transparent  yellow  contour  shows  the  path  of  the  Gulf  Stream  and  is  defined  as 
s/v?  +  v2  >  0.5ms-1.  Strong  negative  fluxes  of  vorticity,  associated  with  a  southward  deflection 
of  the  Gulf  Steam,  are  seen  at  Charleston  Bump  and  after  separation  of  the  current  from  the  coast 
at  Cape  Hatteras.  The  deep  signals  seen  at  Cape  Hatteras  and  along  the  Blake  Ridge  are  due  to  the 
Deep  Western  Boundary  Current  flowing  southward  and  crossing  the  Gulf  Stream  path  right  after 
its  separation. 


13 


Figure  5:  Depth  integrated  eddy  potential  to  eddy  kinetic  energy  conversion  PeKe  =  w'b'  (m2.s~3) 
(left)  and  mean  to  eddy  kinetic  energy  conversion  KmKe  =  HRS  +  VRS  (m2.s~3)  (right)  for  18- 
year  2.5 km  resolution  ROMS  simulation.  Both  energy  conversions  are  at  the  topographic  features, 
especially  the  Charleston  Bump  and  at  the  separation  point  after  Cape  Hatteras. 
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Figure  6:  Comparison  of  lock-exchange  test  solutions  using  explicit,  adaptively  implicit,  and  fully  - 
implicit  vertical  advection  algorithms.  The  setup  is  the  same  as  in  Weak  et  al.(2012),  and  it  is  also 
a  standard  ROMS  test  problem  (Haidvogel  &  Beckmann,  1999)  which  in  its  turn,  is  inspired  by 
the  classical  work  of  Benjamin  (1968).  The  solution  (vertical  along-channel  xz  cross-section  of 
temperature  field)  is  shown  at  17  hours  since  initialization  which  matches  Fig.  2  and  Fig.  5  from 
Weak  et  al.(2012).  The  length  of  the  domain  is  (Hkm,  depth  is  20 m,  grid  resolution  Ax  =  400m, 
Az  =  0.5m.  Explicit  solution  can  be  obtained  for  At  <  60,  beyond  which  the  code  becomes  nu¬ 
merically  unstable.  Note  that  vertical-to-horizontal  aspect  ratios,  both  Az/Ax  =  1/800  -C  1  and 
h/  Ax  =  1/20  <C  1  are  small,  which  means  that  this  grid  does  admit  nonhydrostatic  effects.  Never¬ 
theless,  this  problem  is  known  to  generate  sharp  fronts  with  large  vertical  velocities.  Unlike  Weak 
et  al.(2012)  who  selected  a  Smolarkiewicz  scheme  for  tracer  advection  (the  best  fit  for  this  particu¬ 
lar  problem,  but  is  too  diffusive  for  realistic  long-term  simulations),  we  use  a  third-order  upstream 
scheme  in  the  horizontal,  and  parabolic  spline  in  the  vertical  (for  the  explicit  part)  direction.  They 
also  choose  to  perform  their  tests  with  At  =  Is  resulting  in  vanishingly  small  CFL,  while  our  goal 
here  is  to  push  it  to  the  limit. 
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